# R version 4.0.2 
# packages required

library(dplyr)
library(ggplot2)
library(stargazer)
library(gtsummary)
library(ggpubr)
library(MASS)
library(fauxnaif)
library(vtable)
library(viridis)
library(tmap) 
library(tmaptools)
library(sf)


load("data.RData")
  

######################################################
##################### MAIN PAPER #####################
######################################################


# Descriptive stats of outcome by treatment condition (Table 1)
# Table in paper manually generated based on statistics shown here 

data%>%dplyr::select(`Social distance`, `Policy support`, Trust, Empathy, treatment)%>%filter(treatment!="Control")%>%st(group="treatment")
data%>%dplyr::select(`Social distance`, `Policy support`, Trust, Empathy)%>%st(summ=c("min(x)", "max(x)"))

# Table showing estimated effects of randomized frames (Table 2)

social_distance<-lm(scale(`Social distance`)~treatment, data=subset(data, treatment!="Control"))
policy<-lm(scale(`Policy support`)~treatment, data=subset(data, treatment!="Control"))
empathy<-lm(scale(Empathy)~treatment, data=subset(data, treatment!="Control"))
trust<-lm(scale(Trust)~treatment, data=subset(data, treatment!="Control"))

stargazer(social_distance, policy, trust, empathy, dep.var.labels=c("Social distance", "Policy support", "Trust", "Empathy"), covariate.labels=c("Shared victimization", "Shared goal"), digits=2, omit.stat=c("rsq", "f", "adj.rsq", "ser"), out="main effects.html")

# Graph showing marginal effects of randomized frames for all four outcomes (Figure 1)

matrix_social_distance<-summary(lm(scale(`Social distance`)~treatment, data=subset(data, treatment!="Control")))$coefficients
matrix_social_distance<-matrix_social_distance%>%cbind(c("Social distance", "Social distance", "Social distance"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))
matrix_empathy<-summary(lm(scale(Empathy)~treatment, data=subset(data, treatment!="Control")))$coefficients
matrix_empathy<-matrix_empathy%>%cbind(c("Empathy", "Empathy", "Empathy"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))
matrix_trust<-summary(lm(scale(Trust)~treatment, data=subset(data, treatment!="Control")))$coefficients
matrix_trust<-matrix_trust%>%cbind(c("Trust", "Trust", "Trust"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))
matrix_policy<-summary(lm(scale(`Policy support`)~treatment, data=subset(data, treatment!="Control")))$coefficients
matrix_policy<-matrix_policy%>%cbind(c("Policy support", "Policy support", "Policy support"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))
main_matrix<-rbind(matrix_social_distance, matrix_empathy, matrix_trust, matrix_policy)
main<-as.data.frame(main_matrix, row.names=NULL)

main$V5<-as.factor(main$V5)
main$V5<-factor(main$V5, levels=c("Social distance", "Policy support", "Trust", "Empathy"))
main$`Std. Error`<-as.numeric(main$`Std. Error`)
main$Estimate<-as.numeric(main$Estimate)

limits = aes(ymax = Estimate + (1.96*`Std. Error`), ymin=Estimate - (1.96*`Std. Error`))
dodge = position_dodge(width=0.9)

main%>%filter(V6!="Intercept")%>%ggplot(aes(x = V6, y = Estimate))+facet_wrap(~V5, scales = "free_y", ncol=2)+geom_point(shape=18, size=3)+
  geom_errorbar(limits, position=dodge, width=0.1)+xlab("")+ylab("Change in standard deviation")+theme_classic()+theme(text = element_text(size = 10))+ylim(-0.5,0.5)+coord_flip()+geom_hline(yintercept = 0, color="grey") 


# Subgroup analysis 
# Graph showing the effects on trust and policy support conditional on household income (Figure 2)

matrix_trust_inc<-summary(lm(scale(Trust)~treatment, data=subset(data, treatment!="Control"&income=="Income did not decrease")))$coefficients
matrix_trust_inc<-matrix_trust_inc%>%cbind(c("Income did not decrease", "Income did not decrease", "Income did not decrease"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))

matrix_trust_incde<-summary(lm(scale(Trust)~treatment, data=subset(data, treatment!="Control"&income=="Income decreased")))$coefficients
matrix_trust_incde<-matrix_trust_incde%>%cbind(c("Income decreased", "Income decreased", "Income decreased"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))

trust_inc_matrix<-rbind(matrix_trust_inc,matrix_trust_incde)
trust_inc<-as.data.frame(trust_inc_matrix, row.names=NULL)

trust_inc$V5<-as.factor(trust_inc$V5)
trust_inc$V5<-factor(trust_inc$V5, levels=c("Income did not decrease", "Income decreased"))
trust_inc$`Std. Error`<-as.numeric(trust_inc$`Std. Error`)
trust_inc$Estimate<-as.numeric(trust_inc$Estimate)

matrix_policy_inc<-summary(lm(scale(`Policy support`)~treatment, data=subset(data, treatment!="Control"&income=="Income did not decrease")))$coefficients
matrix_policy_inc<-matrix_policy_inc%>%cbind(c("Income did not decrease", "Income did not decrease", "Income did not decrease"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))

matrix_policy_incde<-summary(lm(scale(`Policy support`)~treatment, data=subset(data, treatment!="Control"&income=="Income decreased")))$coefficients
matrix_policy_incde<-matrix_policy_incde%>%cbind(c("Income decreased", "Income decreased", "Income decreased"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))

policy_inc_matrix<-rbind(matrix_policy_incde,matrix_policy_inc)
policy_inc<-as.data.frame(policy_inc_matrix, row.names=NULL)

policy_inc$V5<-as.factor(policy_inc$V5)
policy_inc$V5<-factor(policy_inc$V5, levels=c("Income did not decrease", "Income decreased"))
policy_inc$`Std. Error`<-as.numeric(policy_inc$`Std. Error`)
policy_inc$Estimate<-as.numeric(policy_inc$Estimate)

trust_inc<-trust_inc%>%cbind(indicator=c("Trust", "Trust", "Trust","Trust", "Trust", "Trust"))
policy_inc<-policy_inc%>%cbind(indicator=c("Policy support", "Policy support", "Policy support","Policy support", "Policy support", "Policy support"))
trust_policy_inc<-rbind(trust_inc, policy_inc)
limits = aes(ymax = Estimate + (1.96*`Std. Error`), ymin=Estimate - (1.96*`Std. Error`))
dodge = position_dodge(width=0.9)

trust_policy_inc%>%filter(V6!="Intercept")%>%
  ggplot(aes(x = V6, y = Estimate, color=V5))+facet_wrap(~indicator, scales = "free_y", ncol=2)+geom_point(shape=18, size=3, position=position_dodge(width=.5))+
  geom_errorbar(limits, position=position_dodge(width=.5), width=0.1)+xlab("")+ylab("Change in standard deviation")+theme_classic()+theme(text = element_text(size = 10))+ylim(-0.5,0.5)+coord_flip()+geom_hline(yintercept = 0, color="grey")+ggtitle('')+guides(color = guide_legend(title = "Subgroup"))+scale_colour_grey()  

# Graph showing the effects on trust and policy support conditional on ethnic majority/minority status (Figure 3)

matrix_trust_bamar<-summary(lm(scale(Trust)~treatment, data=subset(data, treatment!="Control"&bamar=="Bamar")))$coefficients
matrix_trust_bamar<-matrix_trust_bamar%>%cbind(c("Bamar", "Bamar", "Bamar"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))

matrix_trust_em<-summary(lm(scale(Trust)~treatment, data=subset(data, treatment!="Control"&bamar=="Ethnic minorities")))$coefficients
matrix_trust_em<-matrix_trust_em%>%cbind(c("Ethnic minorities", "Ethnic minorities", "Ethnic minorities"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))

trust_bamar_matrix<-rbind(matrix_trust_bamar,matrix_trust_em)
trust_bamar<-as.data.frame(trust_bamar_matrix, row.names=NULL)

trust_bamar$V5<-as.factor(trust_bamar$V5)
trust_bamar$V5<-factor(trust_bamar$V5, levels=c("Ethnic minorities", "Bamar"))
trust_bamar$`Std. Error`<-as.numeric(trust_bamar$`Std. Error`)
trust_bamar$Estimate<-as.numeric(trust_bamar$Estimate)

matrix_policy_bamar<-summary(lm(scale(`Policy support`)~treatment, data=subset(data, treatment!="Control"&bamar=="Bamar")))$coefficients
matrix_policy_bamar<-matrix_policy_bamar%>%cbind(c("Bamar", "Bamar", "Bamar"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))

matrix_policy_em<-summary(lm(scale(`Policy support`)~treatment, data=subset(data, treatment!="Control"&bamar=="Ethnic minorities")))$coefficients
matrix_policy_em<-matrix_policy_em%>%cbind(c("Ethnic minorities", "Ethnic minorities", "Ethnic minorities"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))

policy_bamar_matrix<-rbind(matrix_policy_em,matrix_policy_bamar)
policy_bamar<-as.data.frame(policy_bamar_matrix, row.names=NULL)

policy_bamar$V5<-as.factor(policy_bamar$V5)
policy_bamar$V5<-factor(policy_bamar$V5, levels=c("Ethnic minorities", "Bamar"))
policy_bamar$`Std. Error`<-as.numeric(policy_bamar$`Std. Error`)
policy_bamar$Estimate<-as.numeric(policy_bamar$Estimate)

trust_bamar<-trust_bamar%>%cbind(indicator=c("Trust", "Trust", "Trust","Trust", "Trust", "Trust"))
policy_bamar<-policy_bamar%>%cbind(indicator=c("Policy support", "Policy support", "Policy support","Policy support", "Policy support", "Policy support"))
trust_policy_bamar<-rbind(trust_bamar, policy_bamar)
limits = aes(ymax = Estimate + (1.96*`Std. Error`), ymin=Estimate - (1.96*`Std. Error`))
dodge = position_dodge(width=0.9)

trust_policy_bamar%>%filter(V6!="Intercept")%>%
  ggplot(aes(x = V6, y = Estimate, color=V5))+facet_wrap(~indicator, scales = "free_y", ncol=2)+geom_point(shape=18, size=3, position=position_dodge(width=.5))+
  geom_errorbar(limits, position=position_dodge(width=.5), width=0.1)+xlab("")+ylab("Change in standard deviation")+theme_classic()+theme(text = element_text(size = 10))+ylim(-0.5,0.5)+coord_flip()+geom_hline(yintercept = 0, color="grey")+ggtitle('')+guides(color = guide_legend(title = "Subgroup"))+scale_colour_grey()  


######################################################
##################### APPENDIX B #####################
######################################################

# Geographical distribution of our respondents (Appendix B) 

# Map showing distribution of survey respondents (Figure A2)

map_filt<-data%>%filter(treatment_received==1) 

data_for_map<-map_filt%>%dplyr::select(Tsp_Pcode)%>%filter(!is.na(Tsp_Pcode))%>%group_by(Tsp_Pcode)%>%mutate(count=n())%>%unique()

data_for_map<-data_for_map%>%rename("TS_PCODE"="Tsp_Pcode")

map <- sf::st_read("mmr_polbnda_adm3_mimu_250k.shp", stringsAsFactors = F)
map_and_data<- left_join(map, data_for_map)

map_filtered<-map_and_data%>%
  ggplot() + geom_sf(aes(geometry=geometry), lwd = 0)+geom_point(mapping=aes(geometry=geometry, alpha=0.5,color=count,size=count), stat="sf_coordinates")+ scale_color_continuous(name="", type="viridis", limits=c(1, 260))+ scale_size_continuous(name="", range=c(1,10), limits=c(1,260))+theme(panel.background = element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(),axis.text.x=element_blank(), axis.ticks.y=element_blank(), axis.text.y=element_blank()) + labs(title="No. of respondents from each township")+guides(size="none", alpha="none")

# Distribution of Bamars

map_b<-data%>%filter(treatment_received==1)%>%filter(bamar=="Bamar")

data_b<-map_b%>%dplyr::select(Tsp_Pcode, bamar)%>%filter(!is.na(Tsp_Pcode))%>%group_by(Tsp_Pcode, bamar)%>%mutate(count=n())%>%unique()

data_b<-data_b%>%rename("TS_PCODE"="Tsp_Pcode")

map_and_data_b<- left_join(map, data_b)

map_b<-map_and_data_b%>%
  ggplot() + geom_sf(aes(geometry=geometry), lwd = 0)+geom_point(mapping=aes(geometry=geometry, alpha=0.5,color=count,size=count), stat="sf_coordinates")+ scale_color_continuous(name="", type="viridis", limits=c(1, 260))+ scale_size_continuous(name="", range=c(1,10), limits=c(1,260))+theme(panel.background = element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(),axis.text.x=element_blank(), axis.ticks.y=element_blank(), axis.text.y=element_blank()) + labs(title="No. of Bamar respondents \n from each township")+guides(size="none", alpha="none")

# Distribution of local violence

map_v<-data%>%filter(treatment_received==1)%>%filter(local_violence=="Local violence")

data_v<-map_v%>%dplyr::select(Tsp_Pcode, local_violence)%>%group_by(Tsp_Pcode, local_violence)%>%mutate(count=n())%>%unique()

data_v<-data_v%>%rename("TS_PCODE"="Tsp_Pcode")

map_and_data_v<- left_join(map, data_v)

map_violence<-map_and_data_v%>%
  ggplot() + geom_sf(aes(geometry=geometry), lwd = 0)+geom_point(mapping=aes(geometry=geometry, alpha=0.5, color=count,size=count), stat="sf_coordinates")+ scale_color_continuous(name="", type="viridis", limits=c(1, 260))+ scale_size_continuous(name="", range=c(1,10), limits=c(1,260))+theme(panel.background = element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(),axis.text.x=element_blank(), axis.ticks.y=element_blank(), axis.text.y=element_blank()) + labs(title="Townships with local \n violence (our sample)")+guides(size="none", alpha="none", color="none")

# Distribution of respondents with decreased income

map_inc<-data%>%filter(treatment_received==1)%>%filter(income=="Income decreased")

data_inc<-map_inc%>%dplyr::select(Tsp_Pcode, Township_Name_MMR, Township_Name_Eng, SR_Name_Eng, SR_Pcode, District.SAZ_Pcode, District.SAZ_Name_Eng, income)%>%filter(!is.na(Tsp_Pcode))%>%group_by(Tsp_Pcode, income)%>%mutate(count=n())%>%unique()

colnames(data_inc)[colnames(data_inc)=="Tsp_Pcode"]<-"TS_PCODE"

map_and_data_inc<- left_join(map, data_inc)

map_inc<-map_and_data_inc%>%
  ggplot() + geom_sf(aes(geometry=geometry), lwd = 0)+geom_point(mapping=aes(geometry=geometry, alpha=0.5,color=count,size=count), stat="sf_coordinates")+ scale_color_continuous(name="", type="viridis", limits=c(1, 260))+ scale_size_continuous(name="", range=c(1,10), limits=c(1,260))+theme(panel.background = element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(),axis.text.x=element_blank(), axis.ticks.y=element_blank(), axis.text.y=element_blank()) + labs(title="Respondents with decreased  \n income from each township")+guides(size="none", alpha="none",color="none")

# Maps showing distribution of survey respondents according to ethnicity, local violence and decreased income (Figure A3)

ggarrange(map_b, map_violence, map_inc, nrow=1)


######################################################
##################### APPENDIX C #####################
######################################################

# Our sample versus the national population (Table A1 in Appendix C)
data%>%dplyr::select(age_15to29, urban, Buddhist, states, bamar_num)%>%tbl_summary(missing = "no", statistic=all_continuous()~c("{mean}"), label=list(age_15to29~"Age 15-29 (age 15+)", states~"Reside in state", bamar_num~"Bamar")) # table in paper manually generated based on statistics shown here 

# sample breakdown for less than high school education (>age 25)
data%>%filter(Age_original>3&!is.na(education_num))%>%dplyr::select(less_hs)%>%tbl_summary(missing = "no", statistic=all_continuous()~c("{mean}"), label=less_hs~"Less than high school education (> age 25)") # table in paper manually generated based on statistics shown here 


######################################################
##################### APPENDIX D #####################
######################################################

# Distribution of survey responses by outcomes (Figure A4)

trust_graph<-data%>%filter(!is.na(Trust))%>%group_by(Trust)%>%mutate(count=n())%>%ungroup()%>%mutate(percentage=count/n()*100)%>%group_by(Trust, percentage)%>%summarize%>%
  ggplot(aes(factor(Trust), percentage, label=(paste0(round(percentage, 1), "%"))))+geom_col(position="dodge")+
  geom_text(size=3.5, color="black", position=position_dodge(width=1), vjust=-0.5)+theme_bw()+theme(legend.title = element_blank(), axis.title.x = element_blank())+
  scale_fill_grey(start = 0.3, end = .9)+scale_x_discrete(labels=c("0"="Very unlikely", "1"="Somewhat unlikely", "2"="Somewhat likely", "3"="Very likely"))+
  ggtitle("How likely do you suppose they will return your wallet?")


social_distance_graph<-data%>%filter(!is.na(`Social distance`))%>%group_by(`Social distance`)%>%mutate(count=n())%>%ungroup()%>%mutate(percentage=count/n()*100)%>%group_by(`Social distance`, percentage)%>%summarize%>%
  ggplot(aes(factor(`Social distance`), percentage, label=(paste0(round(percentage, 1), "%"))))+geom_col(position="dodge")+
  geom_text(size=3.5, color="black", position=position_dodge(width=1), vjust=-0.5)+theme_bw()+theme(legend.title = element_blank(), axis.title.x = element_blank())+
  scale_fill_grey(start = 0.3, end = .9)+scale_x_discrete(labels=c("0"="Did not select Rohingya", "1"="Selected Rohingya"))+ggtitle("Which of the following would you not want as a neighbour?")


empathy_graph<-data%>%filter(!is.na(Empathy))%>%group_by(Empathy)%>%mutate(count=n())%>%ungroup()%>%mutate(percentage=count/n()*100)%>%group_by(Empathy, percentage)%>%summarize%>%
  ggplot(aes(factor(Empathy), percentage, label=(paste0(round(percentage, 1), "%"))))+geom_col(position="dodge")+
  geom_text(size=3.5, color="black", position=position_dodge(width=1), vjust=-0.5)+theme_bw()+theme(legend.title = element_blank(), axis.title.x = element_blank())+
  scale_fill_grey(start = 0.3, end = .9)+scale_x_discrete(labels=c("0"="Strongly disagree", "1"="Somewhat disagree", "2"="Somewhat agree", "3"="Strongly agree"))+
  ggtitle("Able to empathize with Rohingya people's sufferings")


policy_graph<-data%>%filter(!is.na(Policy_manual))%>%group_by(Policy_manual)%>%mutate(count=n())%>%ungroup()%>%mutate(percentage=count/n()*100)%>%group_by(Policy_manual, percentage)%>%summarize%>%
  ggplot(aes(factor(Policy_manual), percentage, label=(paste0(round(percentage, 1), "%"))))+geom_col(position="dodge")+
  geom_text(size=3.5, color="black", position=position_dodge(width=1), vjust=-0.5)+theme_bw()+theme(legend.title = element_blank(), axis.title.x = element_blank())+
  scale_fill_grey(start = 0.3, end = .9)+scale_x_discrete(labels=c("1"="No citizenship rights", "2"="Some c rights", "3"="Full c rights", "4"="Other"))+ggtitle("Rohingya people should have")

ggarrange(trust_graph, social_distance_graph, empathy_graph, policy_graph,  nrow=2, ncol=2)

# Item non-response by treatment group (Table A2)

# Table manually generated from the statistics below
# Non-response rate for social distance 
data<-data%>%mutate(social_distance_missing=ifelse((Q57_1=="-99"&Q57_2=="-99"&Q57_3=="-99"&Q57_4=="-99"&Q57_5=="-99"&Q57_6=="-99"&Q57_7=="-99"&Q57_8=="-99"&Q57_9=="-99"&Q57_10=="-99"), 1, ifelse((is.na(Q57_1)|is.na(Q57_2)|is.na(Q57_3)|is.na(Q57_4)|is.na(Q57_5)|is.na(Q57_6)|is.na(Q57_7)|is.na(Q57_8)|is.na(Q57_9)|is.na(Q57_10)), NA, 0)))

data%>%filter(treatment!="Control"&!is.na(social_distance_missing))%>%group_by(treatment, social_distance_missing)%>%summarise(n=n())%>%mutate(freq=n/sum(n)) 

# Non-response rate for policy support
data%>%filter(treatment!="Control"&!is.na(Q66))%>%group_by(treatment, Q66)%>%summarise(n=n())%>%mutate(freq=n/sum(n))%>%filter(Q66=="-99") 

# Non-response rate for trust
data%>%filter(treatment!="Control"&!is.na(Q67))%>%group_by(treatment, Q67)%>%summarise(n=n())%>%mutate(freq=n/sum(n))%>%filter(Q67=="-99") 

# Non-response rate for empathy
data%>%filter(treatment!="Control"&!is.na(Q71))%>%group_by(treatment, Q71)%>%summarise(n=n())%>%mutate(freq=n/sum(n))%>%filter(Q71=="-99") 

# Balance check (Table A3 in Appendix D)

# Table manually generated from the statistics below
data%>%dplyr::select(male, Age_original, urban, education_num, income_level, states, national_id, ethnic_id, Rakhine, bamar_num, local_violence_num, treatment, income_decrease, youth, Buddhist)%>%filter(treatment!="Control")%>%st(group="treatment")

data%>%dplyr::select(male, Age_original, urban, education_num, income_level, states, national_id, ethnic_id, Rakhine, bamar_num, local_violence_num, Buddhist, treatment, income_decrease, youth)%>%filter(treatment!="Control")%>%st(summ=c("min(x)", "max(x)"))

variables_balance<-c("male", "Age_original", "urban", "education_num", "income_level", "states", "national_id", "ethnic_id", "Rakhine", "bamar_num", "local_violence_num", "income_decrease", "youth", "Buddhist")

lm_results <- lapply(variables_balance, function(col){
  lm_formula <- as.formula(paste0(col, "~treatment"))
  lm(lm_formula, data = subset(data, treatment!="Control"))
})

lm_results <- setNames(lm_results, variables_balance) 

lapply(lm_results, function (x) summary(x))

######################################################
##################### APPENDIX E #####################
######################################################

# Assessing treatment delivery (Table A4)

actual_violence_reg<-lm(actual_violence~treatment, data=data)
perceived_violence_reg<-lm(perceived_violence~treatment, data=data)

stargazer(perceived_violence_reg, actual_violence_reg, dep.var.labels=c("Perceived violence", "Violence in Feb 2022"), covariate.labels=c("Placebo", "Shared victimization", "Shared goal"), digits=2, omit.stat=c("rsq", "f", "adj.rsq", "ser"), out="treatment effects.html")

######################################################
##################### APPENDIX F #####################
######################################################

# Alternative reference category (Table A5 in Appendix F)

social_distance_alt<-lm(scale(`Social distance`)~treatment, data=data)
policy_alt<-lm(scale(`Policy support`)~treatment, data=data)
empathy_alt<-lm(scale(Empathy)~treatment, data=data)
trust_alt<-lm(scale(Trust)~treatment, data=data)

stargazer(social_distance_alt, policy_alt, trust_alt, empathy_alt, dep.var.labels=c("Social distance", "Policy support", "Trust", "Empathy"), covariate.labels=c("Placebo", "Shared victimization", "Shared goal"), digits=2, omit.stat=c("rsq", "f", "adj.rsq", "ser"), out="main effects (control as ref).html")


######################################################
##################### APPENDIX G #####################
######################################################

# Robustness check with covariate adjustments (Table A6)

social_distance_covariate<-lm(scale(`Social distance`)~treatment+income_level_imput+states, data=subset(data, treatment!="Control"))
policy_covariate<-lm(scale(`Policy support`)~treatment+income_level_imput+states, data=subset(data, treatment!="Control"))
empathy_covariate<-lm(scale(Empathy)~treatment+income_level_imput+states, data=subset(data, treatment!="Control"))
trust_covariate<-lm(scale(Trust)~treatment+income_level_imput+states, data=subset(data, treatment!="Control"))

stargazer(social_distance_covariate, policy_covariate, trust_covariate, empathy_covariate, dep.var.labels=c("Social distance", "Policy support", "Trust", "Empathy"), covariate.labels=c("Shared victimization", "Shared goal", "Income", "Resides in state"), digits=2, omit.stat=c("rsq", "f", "adj.rsq", "ser"), out="main effects (covariates).html")

# Logistic regression (Table A7)

social_distance_logit<-glm(`Social distance`~relevel(treatment, ref="Placebo"), data=subset(data, treatment!="Control"), family="binomial")
policy_logit<-polr(as.factor(`Policy support`)~relevel(treatment, ref="Placebo"), data=subset(data, treatment!="Control"), Hess=TRUE, method = "logistic")
empathy_logit<-polr(as.factor(Empathy)~relevel(treatment, ref="Placebo"), data=subset(data, treatment!="Control"), Hess=TRUE, method = "logistic")
trust_logit<-polr(as.factor(Trust)~relevel(treatment, ref="Placebo"), data=subset(data, treatment!="Control"), Hess=TRUE, method = "logistic")

stargazer(social_distance_logit, policy_logit, trust_logit, empathy_logit, dep.var.labels=c("Social distance", "Policy support", "Trust", "Empathy"), covariate.labels=c("Shared victimization", "Shared goal"), digits=2, omit.stat=c("rsq", "f", "adj.rsq", "ser", "ll", "aic"), model.names = FALSE, out="main effects (logit).html")

# Alternative coding of policy support (Table A8) 

policy_alt<-lm(scale(as.numeric(policy_alt))~relevel(treatment, ref="Placebo"), data=subset(data, treatment!="Control"))
policy_alt_covariate<-lm(scale(as.numeric(policy_alt))~relevel(treatment, ref="Placebo")+income_level_imput+states, data=subset(data, treatment!="Control"))
policy_alt_logit<-polr(as.factor(policy_alt)~relevel(treatment, ref="Placebo"), data=subset(data, treatment!="Control"), Hess=TRUE, method = "logistic")
policy_alt_logit_covariate<-polr(as.factor(policy_alt)~relevel(treatment, ref="Placebo")+income_level_imput+states, data=subset(data, treatment!="Control"), Hess=TRUE, method = "logistic")

stargazer(policy_alt, policy_alt_covariate, policy_alt_logit, policy_alt_logit_covariate, covariate.labels=c("Shared victimization", "Shared goal", "Income", "Resides in state"), digits=2, omit.stat=c("rsq", "f", "adj.rsq", "ser", "ll", "aic"), dep.var.labels.include = FALSE, out="Policy (alternative coding).html")


######################################################
##################### APPENDIX H #####################
######################################################


# Table showing the effects on trust and policy support conditional on household income (Table A9)

policy_income<-lm(scale(`Policy support`)~relevel(treatment, ref="Placebo")*relevel(income, ref="Income did not decrease"), data=subset(data, treatment!="Control"))
policy_income_logit<-polr(as.factor(`Policy support`)~relevel(treatment, ref="Placebo")*relevel(income, ref="Income did not decrease"), data=subset(data, treatment!="Control"), Hess=TRUE, method = "logistic")
trust_income<-lm(scale(Trust)~relevel(treatment, ref="Placebo")*relevel(income, ref="Income did not decrease"), data=subset(data, treatment!="Control"))
trust_income_logit<-polr(as.factor(Trust)~relevel(treatment, ref="Placebo")*relevel(income, ref="Income did not decrease"), data=subset(data, treatment!="Control"), Hess=TRUE, method = "logistic")


stargazer(policy_income, policy_income_logit, trust_income, trust_income_logit, dep.var.labels.include = FALSE, covariate.labels=c("Shared victimization", "Shared goal", "Decreased income", "Shared victimization x Decreased income", "Shared goal x Decreased income"), digits=2, omit.stat=c("rsq", "f", "adj.rsq", "ser"), column.labels   = c("Policy support", "Trust"),
          column.separate = c(2, 2), model.names = FALSE, out="conditional effects (income).html")


# Table showing the effects on trust and policy support conditional on ethnic majority/minority status (Table A10) 

policy_bamar<-lm(scale(`Policy support`)~relevel(treatment, ref="Placebo")*relevel(bamar, ref="Ethnic minorities"), data=subset(data, treatment!="Control"))
policy_bamar_logit<-polr(as.factor(`Policy support`)~relevel(treatment, ref="Placebo")*relevel(bamar, ref="Ethnic minorities"), data=subset(data, treatment!="Control"), Hess=TRUE, method = "logistic")
trust_bamar<-lm(scale(Trust)~relevel(treatment, ref="Placebo")*relevel(bamar, ref="Ethnic minorities"), data=subset(data, treatment!="Control"))
trust_bamar_logit<-polr(as.factor(Trust)~relevel(treatment, ref="Placebo")*relevel(bamar, ref="Ethnic minorities"), data=subset(data, treatment!="Control"), Hess=TRUE, method = "logistic")


stargazer(policy_bamar, policy_bamar_logit, trust_bamar, trust_bamar_logit, dep.var.labels.include = FALSE, covariate.labels=c("Shared victimization", "Shared goal", "Bamar", "Shared victimization x Bamar", "Shared goal x Bamar"), digits=2, omit.stat=c("rsq", "f", "adj.rsq", "ser"), column.labels   = c("Policy support", "Trust"),
          column.separate = c(2, 2), model.names = FALSE, out="conditional effects (bamar).html")

######################################################
##################### APPENDIX I #####################
######################################################


# Effects on perceived commonality (Figure A5)

matrix_pc<-summary(lm(scale(`Perceived commonalities`)~treatment, data=subset(data, treatment!="Control")))$coefficients
matrix_pc<-matrix_pc%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))
pc_main<-as.data.frame(matrix_pc, row.names=NULL)

pc_main$V5<-as.factor(pc_main$V5)
pc_main$`Std. Error`<-as.numeric(pc_main$`Std. Error`)
pc_main$Estimate<-as.numeric(pc_main$Estimate)

limits = aes(ymax = Estimate + (1.96*`Std. Error`), ymin=Estimate - (1.96*`Std. Error`))
dodge = position_dodge(width=0.9)

pc_main_graph<-pc_main%>%filter(V5!="Intercept")%>%
  ggplot(aes(x = V5, y = Estimate))+geom_point(shape=18, size=3)+
  geom_errorbar(limits, position=dodge, width=0.1)+xlab("")+ylab("Change in standard deviation")+theme_classic()+theme(text = element_text(size = 10))+ylim(-0.5,0.5)+coord_flip()+
  geom_hline(yintercept = 0, color="grey")+ggtitle("Perceived commonalities") 

matrix_pc_inc<-summary(lm(scale(`Perceived commonalities`)~treatment, data=subset(data, treatment!="Control"&income=="Income did not decrease")))$coefficients
matrix_pc_inc<-matrix_pc_inc%>%cbind(c("Income did not decrease", "Income did not decrease", "Income did not decrease"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))

matrix_pc_incde<-summary(lm(scale(`Perceived commonalities`)~treatment, data=subset(data, treatment!="Control"&income=="Income decreased")))$coefficients
matrix_pc_incde<-matrix_pc_incde%>%cbind(c("Income decreased", "Income decreased", "Income decreased"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))

pc_inc_matrix<-rbind(matrix_pc_inc,matrix_pc_incde)
pc_inc<-as.data.frame(pc_inc_matrix, row.names=NULL)

pc_inc$V5<-as.factor(pc_inc$V5)
pc_inc$V5<-factor(pc_inc$V5, levels=c("Income did not decrease", "Income decreased"))
pc_inc$`Std. Error`<-as.numeric(pc_inc$`Std. Error`)
pc_inc$Estimate<-as.numeric(pc_inc$Estimate)

limits = aes(ymax = Estimate + (1.96*`Std. Error`), ymin=Estimate - (1.96*`Std. Error`))
dodge = position_dodge(width=0.9)

pc_inc_graph<-pc_inc%>%filter(V6!="Intercept")%>%
  ggplot(aes(x = V6, y = Estimate, color=V5))+geom_point(shape=18, size=3, position=position_dodge(width=.5))+
  geom_errorbar(limits, position=position_dodge(width=.5), width=0.1)+xlab("")+ylab("")+theme_classic()+theme(text = element_text(size = 10), )+ylim(-0.5,0.5)+coord_flip()+
  geom_hline(yintercept = 0, color="grey")+scale_colour_grey(labels=c("Income did not \ndecrease", "Income decreased"))+guides(color = guide_legend(title = "Subgroup"))

matrix_pc_bamar<-summary(lm(scale(`Perceived commonalities`)~treatment, data=subset(data, treatment!="Control"&bamar=="Bamar")))$coefficients
matrix_pc_bamar<-matrix_pc_bamar%>%cbind(c("Bamar", "Bamar", "Bamar"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))

matrix_pc_em<-summary(lm(scale(`Perceived commonalities`)~treatment, data=subset(data, treatment!="Control"&bamar=="Ethnic minorities")))$coefficients
matrix_pc_em<-matrix_pc_em%>%cbind(c("Ethnic minorities", "Ethnic minorities", "Ethnic minorities"))%>%cbind(c("Intercept", "Shared victimization", "Shared goal"))

pc_bamar_matrix<-rbind(matrix_pc_bamar,matrix_pc_em)
pc_bamar<-as.data.frame(pc_bamar_matrix, row.names=NULL)

pc_bamar$V5<-as.factor(pc_bamar$V5)
pc_bamar$V5<-factor(pc_bamar$V5, levels=c("Ethnic minorities", "Bamar"))
pc_bamar$`Std. Error`<-as.numeric(pc_bamar$`Std. Error`)
pc_bamar$Estimate<-as.numeric(pc_bamar$Estimate)

limits = aes(ymax = Estimate + (1.96*`Std. Error`), ymin=Estimate - (1.96*`Std. Error`))
dodge = position_dodge(width=0.9)

pc_bamar_graph<-pc_bamar%>%filter(V6!="Intercept")%>%
  ggplot(aes(x = V6, y = Estimate, color=V5))+geom_point(shape=18, size=3, position=position_dodge(width=.5))+geom_errorbar(limits, position=position_dodge(width=.5), width=0.1)+xlab("")+
  ylab("")+theme_classic()+theme(text = element_text(size = 10), )+ylim(-0.5,0.5)+coord_flip()+geom_hline(yintercept = 0, color="grey")+
  scale_colour_grey(labels=c("Ethnic minorities", "Bamar"))+guides(color = guide_legend(title = "Subgroup"))

ggarrange(pc_main_graph, pc_inc_graph, pc_bamar_graph, nrow=3, ncol=1)

# Effects on perceived commonality (Table A11)

pc_main<-lm(scale(`Perceived commonalities`)~relevel(treatment, ref="Placebo"), data=subset(data, treatment!="Control"))
pc_main_covariates<-lm(scale(`Perceived commonalities`)~relevel(treatment, ref="Placebo")+income_level_imput+states, data=subset(data, treatment!="Control"))
pc_bamar<-lm(scale(`Perceived commonalities`)~relevel(treatment, ref="Placebo")*relevel(bamar, ref="Ethnic minorities"), data=subset(data, treatment!="Control"))
pc_bamar_covariates<-lm(scale(`Perceived commonalities`)~relevel(treatment, ref="Placebo")*relevel(bamar, ref="Ethnic minorities")+income_level_imput+states, data=subset(data, treatment!="Control"))
pc_income<-lm(scale(`Perceived commonalities`)~relevel(treatment, ref="Placebo")*relevel(income, ref="Income did not decrease"), data=subset(data, treatment!="Control"))
pc_income_covariates<-lm(scale(`Perceived commonalities`)~relevel(treatment, ref="Placebo")*relevel(income, ref="Income did not decrease")+income_level_imput+states, data=subset(data, treatment!="Control"))


stargazer(pc_main, pc_main_covariates, pc_income, pc_income_covariates, pc_bamar, pc_bamar_covariates, covariate.labels=c("Shared victimization", "Shared goal", "Income", "Resides in state", "Decreased income", "V x Income decreased", "G x Income decreased", "Bamar", "V x Bamar", "G x Bamar"), dep.var.labels.include = FALSE, digits=2, omit.stat=c("rsq", "f", "adj.rsq", "ser"), model.names = FALSE, out="Perceived commonalities.html") 

######################################################
##################### APPENDIX J #####################
######################################################

# Effects on co-national outgroups (Table A12)

outgroup_trust<-lm(scale(outgroup_trust)~relevel(treatment, ref="Placebo"), data=subset(data, treatment!="Control"))
federalism<-lm(scale(Federalism)~relevel(treatment, ref="Placebo"), data=subset(data, treatment!="Control"))

stargazer(outgroup_trust, federalism, covariate.labels=c("Shared victimization", "Shared goal"), dep.var.labels=c("Outgroup trust (Co-nationals)", "Support for federalism"), digits=2, omit.stat=c("rsq", "f", "adj.rsq", "ser"), model.names = FALSE, out="Co-nationals.html") 


######################################################
##################### APPENDIX K #####################
######################################################

# Effects in less repressive context (Table A13)

policy_violence<-lm(scale(`Policy support`)~relevel(treatment, ref="Placebo"), data=subset(data, treatment!="Control"&actual_violence==0))
trust_violence<-lm(scale(Trust)~relevel(treatment, ref="Placebo"), data=subset(data, treatment!="Control"&actual_violence==0))

stargazer(policy_violence, trust_violence, dep.var.labels.include = FALSE, covariate.labels=c("Shared victimization", "Shared goal"), digits=2, omit.stat=c("rsq", "f", "adj.rsq", "ser"), column.labels   = c("Policy support", "Trust"), model.names = FALSE, out="lessrepressive.html")



